{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3aacf703",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 前置函数\n",
    "\n",
    "\n",
    "import sympy as sp\n",
    "  \n",
    "# 定义Hilbert矩阵\n",
    "def hilbert_matrix(n):\n",
    "    return sp.Matrix([[sp.Rational(1, i+j-1) for j in range(1, n+1)] for i in range(1, n+1)])\n",
    "\n",
    "# 定义向量s\n",
    "def vector_s(n):\n",
    "    return sp.Matrix([sp.Rational((-1)**j) for j in range(0, n)])\n",
    "\n",
    "# 定义向量t\n",
    "def vector_t(n):\n",
    "    t = []\n",
    "    for j in range(0, n):\n",
    "        sum_val = 0\n",
    "        for i in range(1, j+1):\n",
    "            sum_val += sp.Rational((-1)**(i+j) , i)\n",
    "        t.append(sum_val)\n",
    "    return sp.Matrix(t)\n",
    "\n",
    "# 高斯消元函数\n",
    "def gaussian_elimination(matrix):\n",
    "    n = matrix.shape[0]\n",
    "    for i in range(n):\n",
    "        # 选列主元\n",
    "        max_row = max(range(i, n), key=lambda r: abs(matrix[r, i]))\n",
    "        if i != max_row:\n",
    "            matrix.row_swap(i, max_row)\n",
    "        for j in range(i+1, n):\n",
    "            factor = matrix[j, i] / matrix[i, i]\n",
    "            matrix.row_op(j, lambda x, k: x - factor*matrix[i, k])\n",
    "    # 回代求解\n",
    "    solution = sp.zeros(n, 1)\n",
    "    for i in range(n-1, -1, -1):\n",
    "        solution[i] = matrix[i, -1]\n",
    "        for j in range(i+1, n):\n",
    "            solution[i] -= matrix[i, j] * solution[j]\n",
    "        solution[i] /= matrix[i, i]\n",
    "    return solution\n",
    "\n",
    "def solve_hilbert(b):\n",
    "    n = b.shape[0]\n",
    "    # 创建增广矩阵 [H | b]\n",
    "    augmented_matrix = hilbert_matrix(n).row_join(b)\n",
    "    # 调用高斯消元函数\n",
    "    return gaussian_elimination(augmented_matrix)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4c714683",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n\tbeta=[beta_0,...,beta_n]\n",
      "--------------------------------------------------------------------------------\n",
      "1\t[10, -18]\n",
      "2\t[75, -408, 390]\n",
      "3\t[516, -5700, 13620, -8820]\n",
      "4\t[3405, -63480, 273630, -413280, 202230]\n",
      "5\t[21918, -618870, 4161360, -10780560, 11865420, -4665276]\n",
      "6\t[138775, -5526864, 53241300, -207100320, 379964970, -328592880, 107975868]\n"
     ]
    }
   ],
   "source": [
    "# beta的精确求解\n",
    "\n",
    "# 打印表头\n",
    "print(\"n\\tbeta=[beta_0,...,beta_n]\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "for n in range(1, 7):\n",
    "    beta = solve_hilbert(vector_s(n+1))\n",
    "    # 打印 beta\n",
    "    beta = [f\"{val}\" for val in beta]\n",
    "    print(f\"{n}\\t[{', '.join(beta)}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "31430af8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n\tgamma=[gamma_0,...,gamma_n]\n",
      "--------------------------------------------------------------------------------\n",
      "1\t[-6, 12]\n",
      "2\t[-51, 282, -270]\n",
      "3\t[-1070/3, 3950, -9440, 18340/3]\n",
      "4\t[-14155/6, 44000, -189665, 859390/3, -140175]\n",
      "5\t[-75957/5, 428967, -2884434, 7472514, -8224482, 16168614/5]\n",
      "6\t[-192381/2, 19154646/5, -36904056, 143551002, -263371647, 227763228, -374215842/5]\n"
     ]
    }
   ],
   "source": [
    "# gamma的精确求解\n",
    "\n",
    "# 打印表头\n",
    "print(\"n\\tgamma=[gamma_0,...,gamma_n]\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "for n in range(1, 7):\n",
    "    gamma = solve_hilbert(vector_t(n+1))\n",
    "    # 打印 gamma\n",
    "    gamma = [f\"{val}\" for val in gamma]\n",
    "    print(f\"{n}\\t[{', '.join(gamma)}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "56f3ea90",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n\talpha=[alpha_0,...,alpha_n]\n",
      "--------------------------------------------------------------------------------\n",
      "1\t[-6 + 10*log(2), 12 - 18*log(2)]\n",
      "2\t[-51 + 75*log(2), 282 - 408*log(2), -270 + 390*log(2)]\n",
      "3\t[-1070/3 + 516*log(2), 3950 - 5700*log(2), -9440 + 13620*log(2), 18340/3 - 8820*log(2)]\n",
      "4\t[-14155/6 + 3405*log(2), 44000 - 63480*log(2), -189665 + 273630*log(2), 859390/3 - 413280*log(2), -140175 + 202230*log(2)]\n",
      "5\t[-75957/5 + 21918*log(2), 428967 - 618870*log(2), -2884434 + 4161360*log(2), 7472514 - 10780560*log(2), -8224482 + 11865420*log(2), 16168614/5 - 4665276*log(2)]\n",
      "6\t[-192381/2 + 138775*log(2), 19154646/5 - 5526864*log(2), -36904056 + 53241300*log(2), 143551002 - 207100320*log(2), -263371647 + 379964970*log(2), 227763228 - 328592880*log(2), -374215842/5 + 107975868*log(2)]\n"
     ]
    }
   ],
   "source": [
    "# alpha的精确求解\n",
    "\n",
    "# 打印表头\n",
    "print(\"n\\talpha=[alpha_0,...,alpha_n]\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "for n in range(1, 7):\n",
    "    beta = solve_hilbert(vector_s(n+1))\n",
    "    gamma = solve_hilbert(vector_t(n+1))\n",
    "    alpha = beta * sp.log(2) + gamma\n",
    "    # 打印 alpha\n",
    "    alpha = [f\"{val}\" for val in alpha]\n",
    "    print(f\"{n}\\t[{', '.join(alpha)}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "10d9b3e7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n\talpha=[alpha_0,...,alpha_n]\n",
      "--------------------------------------------------------------------------------\n",
      "1\t[0.931471805599453, -0.476649250079015]\n",
      "2\t[0.986038541995896, -0.804049668457651, 0.327400418378659]\n",
      "3\t[0.997278502265090, -0.938929191688203, 0.664599226454811, -0.224799205384443]\n",
      "4\t[0.999483139947188, -0.983021945328801, 0.863016617833637, -0.533448480884545, 0.154324637725949]\n",
      "5\t[0.999903512880337, -0.995633133337833, 0.951294933911413, -0.768857323564589, 0.419159585610032, -0.105933979619294]\n",
      "6\t[0.999982206412824, -0.998938261065632, 0.984346218407154, -0.901062458753586, 0.667044192552567, -0.324072420597076, 0.0727128088474274]\n"
     ]
    }
   ],
   "source": [
    "# alpha机器精度的求解\n",
    "\n",
    "import math\n",
    "\n",
    "# 打印表头\n",
    "print(\"n\\talpha=[alpha_0,...,alpha_n]\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "for n in range(1, 7):\n",
    "    beta = solve_hilbert(vector_s(n+1))\n",
    "    gamma = solve_hilbert(vector_t(n+1))\n",
    "    alpha = beta * math.log(2) + gamma\n",
    "    # 打印 alpha\n",
    "    alpha = [f\"{val}\" for val in alpha]\n",
    "    print(f\"{n}\\t[{', '.join(alpha)}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c241aa44",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n\talpha=[alpha_0,...,alpha_n]\n",
      "--------------------------------------------------------------------------------\n",
      "1\t[0.931500000000000, -0.476700000000001]\n",
      "2\t[0.986250000000004, -0.805200000000020, 0.328500000000019]\n",
      "3\t[0.998733333333353, -0.955000000000209, 0.703000000000491, -0.249666666666982]\n",
      "4\t[1.00908333333342, -1.16200000000147, 1.63450000000618, -1.69866666667584, 0.724500000004428]\n",
      "5\t[1.06170000000045, -2.74050000001264, 12.6840000000844, -31.1640000002177, 33.8730000002389, -13.2594000000938]\n",
      "6\t[1.39125000000332, -16.5816000001341, 151.095000001301, -584.808000005086, 1071.95550000937, -926.772000008126, 304.504200002677]\n"
     ]
    }
   ],
   "source": [
    "# 取ln2为0.69315时, alpha的求解\n",
    "\n",
    "# 打印表头\n",
    "print(\"n\\talpha=[alpha_0,...,alpha_n]\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "for n in range(1, 7):\n",
    "    r = vector_s(n+1) * 0.69315 + vector_t(n+1)\n",
    "    alpha = solve_hilbert(r)\n",
    "    # 打印 alpha\n",
    "    alpha = [f\"{val}\" for val in alpha]\n",
    "    print(f\"{n}\\t[{', '.join(alpha)}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "668e38fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n\talpha=[alpha_0,...,alpha_n]\n",
      "--------------------------------------------------------------------------------\n",
      "1\t[0.931499990579881, -0.476699982477002]\n",
      "2\t[0.986248358872464, -0.805190837280058, 0.328491164987320]\n",
      "3\t[0.998415516878830, -0.951423697781147, 0.694392620130081, -0.244071631181087]\n",
      "4\t[0.994861561853290, -0.893606956770126, 0.471412836652721, 0.0631042699867720, -0.139109121889517]\n",
      "5\t[0.995675974159657, -0.908098652527691, 0.523036319112395, 0.0247748214019527, -0.182562221696320, 0.0449795044140929]\n",
      "6\t[0.997987897852757, -0.938697151663587, 0.582433501229182, 0.0801931619085124, -0.250541338216056, -0.165637605572610, 0.196226561793087]\n"
     ]
    }
   ],
   "source": [
    "# 取ln2为0.69315时, alpha的Tikhonov正则化求解\n",
    "\n",
    "# 打印表头\n",
    "print(\"n\\talpha=[alpha_0,...,alpha_n]\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "#Tikhonov正则化参数\n",
    "Alpha = 1e-10\n",
    "\n",
    "for n in range(1, 7):\n",
    "    H = hilbert_matrix(n+1)\n",
    "    H_trans = H.transpose()\n",
    "    A = Alpha * sp.eye(n+1) + H_trans * H\n",
    "    \n",
    "    r = vector_s(n+1) * 0.69315 + vector_t(n+1)\n",
    "    b = H_trans * r\n",
    "    # 创建增广矩阵 [A | b]\n",
    "    augmented_matrix = A.row_join(b)\n",
    "    # 调用高斯消元函数\n",
    "    alpha = gaussian_elimination(augmented_matrix)\n",
    "    # 打印 alpha\n",
    "    alpha = [f\"{val}\" for val in alpha]\n",
    "    print(f\"{n}\\t[{', '.join(alpha)}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44ff60cf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
